%%
%Example waveforms and FRA for Figure 3D
load('example_FRA')
figure
imagesc(example_FRA)
colormap('hot')
ax = gca;
ax.YTick = 1:5; ax.YTickLabel = {'80' '65' '50' '35' '20'};
ax.XTick = 1:8; ax.XTickLabel = {'4' '5.7' '8' '11.3' '16' '22.6' '32' '45.3'};
xlabel('Frequency (kHz)')
ylabel('Intensity (dB SPL)')


load('example_calcium_traces')
load('example_spike_traces')
intensities = [20:15:80];
figure(2); figure(3);
for j = 1:5 %Loop through intensities
    
    figure(2)
    subplot(1,5,j)
    plot(calcium_traces(:,j))
    hold on
    plot([15 15],[-50 500],'k--')
    box off
    xlabel('Time (sec)')
    ylabel('Response (dF/F %)')
    ax = gca; ax.XTick = [0 30 60 90]; ax.XTickLabel = {'0' '1' '2' '3'};
    ylim([-50 500])
    title(num2str(intensities(j)))

    figure(3)
    subplot(1,5,j)
    plot(spike_traces(:,j))
    hold on
    plot([15 15],[-1 125],'k--')
    box off
    xlabel('Time (sec)')
    ylabel('Response (events/sec)')
    ax = gca; ax.XTick = [0 30 60 90]; ax.XTickLabel = {'0' '1' '2' '3'};
    ylim([-1 125])
    title(num2str(intensities(j)))

end

%%
%Next make the representative maps in Figure 3E
load('sham_example_mapping.mat')
jet = jet;
colors = jet(1:8:64,:);
days = {'Day -4' 'Day 0' 'Day 21'};
    
for i = 1:3 %Loop through days

    figure;
    h = imshow(mapData(i).im',[]);
    set(h,'AlphaData',0.8)
    hold on

    for n=1:length(mapData(i).BFs) %Loop through cells
        color = colors(mapData(i).BFs(n),:); %Get correct color
        plot(mapData(i).allBounds{n}(:,1),mapData(i).allBounds{n}(:,2),'linewidth',2,'color',color);         
    end

    axis square; title(days{i})
    %Add svm line
    plot(mapData(i).y,mapData(i).x,'--k')
    ylim([0 512])
    xlim([0 512])
end


load('trauma_example_mapping.mat')
jet = jet;
colors = jet(1:8:64,:);
days = {'Day -4' 'Day 0' 'Day 21'};
    
for i = 1:3 %Loop through days

    figure;
    h = imshow(mapData(i).im',[]);
    set(h,'AlphaData',0.8)
    hold on

    for n=1:length(mapData(i).BFs) %Loop through cells
        color = colors(mapData(i).BFs(n),:); %Get correct color
        plot(mapData(i).allBounds{n}(:,1),mapData(i).allBounds{n}(:,2),'linewidth',2,'color',color);         
    end

    axis square; title(days{i})
    %Add svm line
    plot(mapData(i).y,mapData(i).x,'--k')
    ylim([0 512])
    xlim([0 512])
end

%%
%Log scale plot for median BF's, Figure 3G
load('BF_Info')
%Median BF plot
BF_noise = cell(1,4); BF_sham = cell(1,4);
for j = 1:size(BF_Info,2) %Loop through days
    for i = 1:size(BF_Info,1) %Loop through mice
        if ismember(i,[1 3 4 6])
            BF_noise{j} = [BF_noise{j}; BF_Info{i,j}];
        elseif ismember(i,[2 5 7 8])
            BF_sham{j} = [BF_sham{j}; BF_Info{i,j}];
        end
    end
end

colors = {[0 0 0] [0.35 0.35 0.35] [0.7 0.7 0.7]};
threshold = 0.8; dist_thresh = 500;

figure
title('Noise Exposed Mice')
for day = 1:3
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Change the line below  to BF_noise to get the noise-exposed mouse plot,
%etc.
BFs = BF_sham{day}(:,1); dprimes = BF_sham{day}(:,2); dists = BF_sham{day}(:,3);

%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = [];

%Now sort the data
[dists,inds] = sort(dists); BFs = BFs(inds); nGroups = 5;
dists = dists.*1.6328; %To convert to microns

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = [];

%Next go through the neurons and get average distance and response measures
num1 = ceil(length(BFs)/nGroups); num2 = mod(length(BFs), nGroups);
groups = ones(1,nGroups).*num1; 
if num2 ~= 0; groups(1:(nGroups-num2)) = groups(1:(nGroups-num2)) - 1; end

clearvars gBF gDists gErr
start = 0; stop = 0;
for i = 1:nGroups
    start = stop + 1; stop = stop + groups(i);
    gBF(i) = exp(median(log(BFs(start:stop))));
    gDists(i) = mean(dists(start:stop));
    
    %To get the error, use a bootstrap
    median_sims = zeros(1,1000);
    for k = 1:10000
        sample = randsample(BFs(start:stop)',(stop-start+1),true);
        median_sims(k) = exp(median(log(sample)));
    end
    gErr(i) = std(median_sims);
end

%Get logged error
gBF = gBF./1000; gErr = gErr./1000;
gErr_up = log2(gBF+gErr) - log2(gBF);
gErr_down = log2(gBF) - log2(gBF - gErr);
newErr = [gErr_up;gErr_down];

if day == 1
errorbar(gDists,log2(gBF),gErr_up,gErr_down,'Color',colors{day},'LineStyle','--','CapSize',0)
hold on
p{day} = plot(gDists,log2(gBF),'Color',colors{day},'Marker','.','MarkerSize',20,'LineStyle','none');
else
errorbar(gDists,log2(gBF),gErr_up,gErr_down,'Color',colors{day},'CapSize',0)
hold on
p{day} = plot(gDists,log2(gBF),'Color',colors{day},'LineWidth',2,'Marker','.','MarkerSize',20);    
end

end
box off
xlabel('Distance to SVM Line (Microns)')
ylabel('Median BF (kHz)')
legend([p{1} p{2} p{3}], {'Baseline' 'Day 0' 'Day 21'})
ylim([2 5])
xlim([-450 450])
ax = gca;
ax.YTick = [2 3 4 5];
ax.YTickLabel = {'4' '8' '16' '32'};

%%
%Now plot Figure 3H
load('BF_Info')
BF_noise = cell(1,4); BF_sham = cell(1,4);
for j = 1:size(BF_Info,2) %Loop through days
    for i = 1:size(BF_Info,1) %Loop through mice
        if ismember(i,[1 3 4 6])
            BF_noise{j} = [BF_noise{j}; BF_Info{i,j}];
        elseif ismember(i,[2 5 7 8])
            BF_sham{j} = [BF_sham{j}; BF_Info{i,j}];
        end
    end
end


%Decide on parameters
threshold = 0.8; dist_thresh = 500;
lowBound = 9000; highBound = 17000; %Only include 11 and 16 kHz

for day = 1:3
    
BFs = BF_noise{day}(:,1); dprimes = BF_noise{day}(:,2); dists = BF_noise{day}(:,3);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = [];

%Now sort the data
[dists,inds] = sort(dists); BFs = BFs(inds); nGroups = 5;
dists = dists.*1.6328; %To convert to microns

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = [];  

%Now do a calculation for deafferented vs spared. Also need to get a
%measurement of error
for k = 1:2 %k is 1 is spared
    
    if k == 1
        
        inds = dists < 0; %Negative is spared
        percents_noise(k,day) = mean(BFs(inds) > lowBound & BFs(inds) < highBound);
        
        percents_sim = zeros(1,1000);
        for n = 1:1000 %Bootstrap error here
            sample = randsample(BFs(inds),sum(inds),'true');
            percents_sim(n) = mean(sample > lowBound & sample < highBound);
        end
        errors_noise(k,day) = std(percents_sim);
        
    elseif k == 2
        
        inds = dists > 0; %Negative is spared
        percents_noise(k,day) = mean(BFs(inds) > lowBound & BFs(inds) < highBound);
        
        percents_sim = zeros(1,1000);
        for n = 1:1000 %Bootstrap error here
            sample = randsample(BFs(inds),sum(inds),'true');
            percents_sim(n) = mean(sample > lowBound & sample < highBound);
        end
        errors_noise(k,day) = std(percents_sim);
        
        
    end
    
    
    
end


BFs = BF_sham{day}(:,1); dprimes = BF_sham{day}(:,2); dists = BF_sham{day}(:,3);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = [];

%Now sort the data
[dists,inds] = sort(dists); BFs = BFs(inds); nGroups = 5;
dists = dists.*1.6328; %To convert to microns

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = [];  

%Now do a calculation for deafferented vs spared
for k = 1:2 %k is 1 is spared
    
    
    if k == 1
        
        inds = dists < 0; %Negative is spared
        percents_sham(k,day) = mean(BFs(inds) > lowBound & BFs(inds) < highBound);
        
        percents_sim = zeros(1,1000);
        for n = 1:1000 %Bootstrap error here
            sample = randsample(BFs(inds),sum(inds),'true');
            percents_sim(n) = mean(sample > lowBound & sample < highBound);
        end
        errors_sham(k,day) = std(percents_sim);
        
    elseif k == 2
        
        inds = dists > 0; %Negative is spared
        percents_sham(k,day) = mean(BFs(inds) > lowBound & BFs(inds) < highBound);
        
        percents_sim = zeros(1,1000);
        for n = 1:1000 %Bootstrap error here
            sample = randsample(BFs(inds),sum(inds),'true');
            percents_sim(n) = mean(sample > lowBound & sample < highBound);
        end
        errors_sham(k,day) = std(percents_sim);
        
    end
    
    
end



end


figure
subplot(1,2,1)
errorbar(percents_noise(1,:),errors_noise(1,:),'.-','Color','b','LineWidth',2,'MarkerSize',20,'CapSize',0)
hold on
errorbar(percents_sham(1,:),errors_sham(1,:),'.-','Color','k','LineWidth',2,'MarkerSize',20,'CapSize',0)
ylim([0 1])
box off
ax = gca;
ax.XTick = [1:3]; ax.XTickLabel = {'Baseline' 'Day 0' 'Day 21'};
xlabel('Day Re: Exposure')
ylabel('Percent Edge Frequency BF')
title('Intact Region')


subplot(1,2,2)
errorbar(percents_noise(2,:),errors_noise(2,:),'.-','Color','b','LineWidth',2,'MarkerSize',20,'CapSize',0)
hold on
errorbar(percents_sham(2,:),errors_sham(2,:),'.-','Color','k','LineWidth',2,'MarkerSize',20,'CapSize',0)
ylim([0 1])
box off
ax = gca;
ax.XTick = [1:3]; ax.XTickLabel = {'Baseline' 'Day 0' 'Day 21'};
xlabel('Day Re: Exposure')
ylabel('Percent Edge Frequency BF')
title('Deafferented Region')
legend({'Noise' 'Sham'})


%%
%Plot in a similar style to above but now for threshold. Figure 3I
load('BF_Info')
BF_noise = cell(1,4); BF_sham = cell(1,4);
for j = 1:size(BF_Info,2) %Loop through days
    for i = 1:size(BF_Info,1) %Loop through mice
        if ismember(i,[1 3 4 6])
            BF_noise{j} = [BF_noise{j}; BF_Info{i,j}];
        elseif ismember(i,[2 5 7 8])
            BF_sham{j} = [BF_sham{j}; BF_Info{i,j}];
        end
    end
end


threshold = 0.5; dist_thresh = 500;

for day = 1:3
    
BFs = BF_sham{day}(:,1); dprimes = BF_sham{day}(:,2); dists = BF_sham{day}(:,3);
threshs = BF_sham{day}(:,4);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = []; threshs(to_delete) = [];


%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = []; threshs(idxs) = [];

for k = 1:2 %For different regions
    
    if k == 1 %Spared
        idxs = dists < 0; %Negative distances is spared
        threshs_sham(k,day) = mean(threshs(idxs));
        errors_sham(k,day) = std(threshs(idxs))./sqrt(sum(idxs));       
    elseif k ==2 
        idxs = dists > 0;
        threshs_sham(k,day) = mean(threshs(idxs));
        errors_sham(k,day) = std(threshs(idxs))./sqrt(sum(idxs));  
    end
    
end




%Same as above but now for the noise
BFs = BF_noise{day}(:,1); dprimes = BF_noise{day}(:,2); dists = BF_noise{day}(:,3);
threshs = BF_noise{day}(:,4);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = []; threshs(to_delete) = [];

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = []; threshs(idxs) = [];

for k = 1:2 %For different regions
    
    if k == 1 %Spared
        idxs = dists < 0; %Negative distances is spared
        threshs_noise(k,day) = mean(threshs(idxs));
        errors_noise(k,day) = std(threshs(idxs))./sqrt(sum(idxs));       
    elseif k ==2 
        idxs = dists > 0;
        threshs_noise(k,day) = mean(threshs(idxs));
        errors_noise(k,day) = std(threshs(idxs))./sqrt(sum(idxs));  
    end
    
end
        

end



figure
subplot(1,2,1)
errorbar(threshs_noise(1,:),errors_noise(1,:),'.-','Color','b','LineWidth',2,'MarkerSize',20,'CapSize',0)
hold on
errorbar(threshs_sham(1,:),errors_sham(1,:),'.-','Color','k','LineWidth',2,'MarkerSize',20,'CapSize',0)
ylim([20 40])
box off
ax = gca;
ax.XTick = [1:3]; ax.XTickLabel = {'Baseline' 'Day 0' 'Day 21'};
xlabel('Day Re: Exposure')
ylabel('Mean Threshold')
title('Intact Region')


subplot(1,2,2)
errorbar(threshs_noise(2,:),errors_noise(2,:),'.-','Color','b','LineWidth',2,'MarkerSize',20,'CapSize',0)
hold on
errorbar(threshs_sham(2,:),errors_sham(2,:),'.-','Color','k','LineWidth',2,'MarkerSize',20,'CapSize',0)
ylim([20 40])
box off
ax = gca;
ax.XTick = [1:3]; ax.XTickLabel = {'Baseline' 'Day 0' 'Day 21'};
xlabel('Day Re: Exposure')
ylabel('Mean Threshold')
title('Deafferented Region')
legend({'Noise' 'Sham'})

%%
%Finally the statistics for Figure 3H, 3I
load('BF_Info') 

BF_noise = cell(1,4); BF_sham = cell(1,4);
for j = 1:size(BF_Info,2) %Loop through days
    for i = 1:size(BF_Info,1) %Loop through mice
        if ismember(i,[1 3 4 6])
            BF_noise{j} = [BF_noise{j}; BF_Info{i,j}];
        elseif ismember(i,[2 5 7 8])
            BF_sham{j} = [BF_sham{j}; BF_Info{i,j}];
        end
    end
end


%Decide on parameters
threshold = 0.8; dist_thresh = 500;
lowBound = 9000; highBound = 17000; %Only include 11 and 16 kHz
allBFValues = []; treatment = []; region = []; timePoint = [];
%For treatment, 1 is exposure and 2 is sham. For region, 1 is spared and 2
%is deafferented. For time, 1 is the first day, etc.

for day = 1:3
    
BFs = BF_noise{day}(:,1); dprimes = BF_noise{day}(:,2); dists = BF_noise{day}(:,3);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = [];

dists = dists.*1.6328; %To convert to microns

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = [];  

%Now do a calculation for deafferented vs spared. Also need to get a
%measurement of error
for k = 1:2 %k is 1 is spared
    if k == 1
        
        inds = dists < 0; %Negative is spared
        currBFs = BFs(inds); currBFs = double(currBFs > lowBound & currBFs < highBound);
        %Update variables
        allBFValues = [allBFValues; currBFs]; treatment = [treatment; ones(length(currBFs),1)];
        region = [region; ones(length(currBFs),1).*k]; timePoint = [timePoint; ones(length(currBFs),1).*day];
               
    elseif k == 2
        
        inds = dists > 0; %Negative is spared
        currBFs = BFs(inds); currBFs = double(currBFs > lowBound & currBFs < highBound);
        %Update variables
        allBFValues = [allBFValues; currBFs]; treatment = [treatment; ones(length(currBFs),1)];
        region = [region; ones(length(currBFs),1).*k]; timePoint = [timePoint; ones(length(currBFs),1).*day];

    end
end


BFs = BF_sham{day}(:,1); dprimes = BF_sham{day}(:,2); dists = BF_sham{day}(:,3);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = [];

%Now sort the data
dists = dists.*1.6328; %To convert to microns

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = [];  

%Now do a calculation for deafferented vs spared
for k = 1:2 %k is 1 is spared
    if k == 1
        
        inds = dists < 0; %Negative is spared
        currBFs = BFs(inds); currBFs = double(currBFs > lowBound & currBFs < highBound);
        %Update variables
        allBFValues = [allBFValues; currBFs]; treatment = [treatment; ones(length(currBFs),1).*2];
        region = [region; ones(length(currBFs),1).*k]; timePoint = [timePoint; ones(length(currBFs),1).*day];
        
    elseif k == 2
        
        inds = dists > 0; %Negative is spared
        currBFs = BFs(inds); currBFs = double(currBFs > lowBound & currBFs < highBound);
        %Update variables
        allBFValues = [allBFValues; currBFs]; treatment = [treatment; ones(length(currBFs),1).*2];
        region = [region; ones(length(currBFs),1).*k]; timePoint = [timePoint; ones(length(currBFs),1).*day];
        
    end      
end
end

%Now perform the ANOVA
[p,tbl,stats,terms] = anovan(allBFValues,{treatment region timePoint},'model','full',...
    'varnames',{'Treatment' 'Region' 'Time'});

%Also do a version with two separate two-way ANOVA's for each region.
indexes = find(region == 1); %Looking at spared region
allBF_spare = allBFValues(indexes);
treat_spare = treatment(indexes);
time_spare = timePoint(indexes);
[p,tbl,stats,terms] = anovan(allBF_spare,{treat_spare time_spare},'model','full',...
    'varnames',{'Treatment' 'Time'});

indexes = find(region == 2); %Looking at deafferented region
allBF_deaff = allBFValues(indexes);
treat_deaff = treatment(indexes);
time_deaff = timePoint(indexes);
[p,tbl,stats,terms] = anovan(allBF_deaff,{treat_deaff time_deaff},'model','full',...
    'varnames',{'Treatment' 'Time'});


%%
%Now do the same approach but now for threshold.

%Plot in a similar style to above but now for threshold
BF_noise = cell(1,4); BF_sham = cell(1,4);
for j = 1:size(BF_Info,2) %Loop through days
    for i = 1:size(BF_Info,1) %Loop through mice
        if ismember(i,[1 3 4 6])
            BF_noise{j} = [BF_noise{j}; BF_Info{i,j}];
        elseif ismember(i,[2 5 7 8])
            BF_sham{j} = [BF_sham{j}; BF_Info{i,j}];
        end
    end
end

%Contraints
threshold = 0.5; dist_thresh = 500;
allThreshValues = []; treatment = []; region = []; timePoint = [];
%For treatment, 1 is exposure and 2 is sham. For region, 1 is spared and 2
%is deafferented. For time, 1 is the first day, etc.

for day = 1:3
    
BFs = BF_sham{day}(:,1); dprimes = BF_sham{day}(:,2); dists = BF_sham{day}(:,3);
threshs = BF_sham{day}(:,4);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = []; threshs(to_delete) = [];

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = []; threshs(idxs) = [];

for k = 1:2 %For different regions
    
    if k == 1 %Spared
        idxs = dists < 0; %Negative distances is spared
        currThresh = threshs(idxs); 
        %Update variables
        allThreshValues = [allThreshValues; currThresh]; treatment = [treatment; ones(length(currThresh),1)];
        region = [region; ones(length(currThresh),1).*k]; timePoint = [timePoint; ones(length(currThresh),1).*day];
     
    elseif k ==2 
        idxs = dists > 0;
        currThresh = threshs(idxs); 
        %Update variables
        allThreshValues = [allThreshValues; currThresh]; treatment = [treatment; ones(length(currThresh),1)];
        region = [region; ones(length(currThresh),1).*k]; timePoint = [timePoint; ones(length(currThresh),1).*day]; 
    end
    
end




%Same as above but now for the noise
BFs = BF_noise{day}(:,1); dprimes = BF_noise{day}(:,2); dists = BF_noise{day}(:,3);
threshs = BF_noise{day}(:,4);
%Get rid of non well-tuned cells
to_delete = dprimes < threshold;
BFs(to_delete) = []; dists(to_delete) = []; threshs(to_delete) = [];

%Get rid of cells at certain distance away
idxs = abs(dists) > dist_thresh;
dists(idxs) = []; BFs(idxs) = []; threshs(idxs) = [];

for k = 1:2 %For different regions
    
    if k == 1 %Spared
        idxs = dists < 0; %Negative distances is spared
        currThresh = threshs(idxs); 
        %Update variables
        allThreshValues = [allThreshValues; currThresh]; treatment = [treatment; ones(length(currThresh),1).*2];
        region = [region; ones(length(currThresh),1).*k]; timePoint = [timePoint; ones(length(currThresh),1).*day];      
    elseif k ==2 
        idxs = dists > 0;
        currThresh = threshs(idxs); 
        %Update variables
        allThreshValues = [allThreshValues; currThresh]; treatment = [treatment; ones(length(currThresh),1).*2];
        region = [region; ones(length(currThresh),1).*k]; timePoint = [timePoint; ones(length(currThresh),1).*day]; 
    end
    
end
        

end

%Now perform the ANOVA
[p,tbl,stats,terms] = anovan(allThreshValues,{treatment region timePoint},'model','full',...
    'varnames',{'Treatment' 'Region' 'Time'});

%Also do a version with two separate two-way ANOVA's for each region.
indexes = find(region == 1); %Looking at spared region
allThresh_spare = allThreshValues(indexes);
treat_spare = treatment(indexes);
time_spare = timePoint(indexes);
[p,tbl,stats,terms] = anovan(allThresh_spare,{treat_spare time_spare},'model','full',...
    'varnames',{'Treatment' 'Time'});

indexes = find(region == 2); %Looking at deafferented region
allThresh_deaff = allThreshValues(indexes);
treat_deaff = treatment(indexes);
time_deaff = timePoint(indexes);
[p,tbl,stats,terms] = anovan(allThresh_deaff,{treat_deaff time_deaff},'model','full',...
    'varnames',{'Treatment' 'Time'});
